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Multi-particle collision dynamics (MPCD) is a mesoscopic simulation method 
to simulate fluid particle-like flows. MPCD has been widely used to simulate 
various problems in condensed matter. In this study, hydrostatic behavior of 
gas in the Earth’s atmospheric layer is simulated by using MPCD method. 
The simulation is carried out by assuming the system under ideal state and 
is affected only by gravitational force. Gas particles are homogeneous and 
placed in 2D box. Interaction of the particles with the box is applied through 
implementation of boundary conditions (BC). Periodic BC is applied on the 
left and the right side, specular reflection on the top side, while bounce- 
back on the bottom side. Simulation program is executed in Arch Linux and 
running in notebook with processor Intel i5 @2700 MHz with 10 GB DDR3 
RAM. The results show behaviors of the particles obey kinetic theory for 
ideal gas when gravitational acceleration value is proportional to the particle 
mass. Density distribution as a function of altitude also meets atmosphere’s 
hydrostatic theory. 


1 Introduction 

Atmosphere is a layer that surrounds the earth and is mostly composed of molecular gases 
and plasma in upper atmosphere. Atmosphere is interesting to study since it has a lot of 
complex phenomena and still under active study. The complexity is mostly contributed 
by the interaction of fluids, such as the interaction between ocean and neutral gas, or 
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between neutral gas and plasma. Hence, any computational fluid dynamics methods can 
be used to observe this complexity. 

In the other hand, Multi-particle collision dynamics (MPCD) that was first proposed 
by Malevanets and Kapral in 1999 [8] offers a simple computational method to model 
fluids dynamics. In MPCD method, as well as in the other particle-based methods, 
the fluids are assumed as a collection of particle-like objects with zero volume. MPCD 
algorithm is similar to Direct Simulation Monte Carlo (DSMC) [2], the difference is 
only on the collision mechanism among the particles. The collision in DSMC method is 
represented by binary-collision between the pair of particles; while in MPCD, the collision 
is represented by stochastic rotation dynamics. MPCD with this collision mechanism 
is known as Stochastic Rotation Dynamics (SRD). MPCD with Anderson Thermostat 
(MPCD-AT) can be used as an alternative to SRD mechanism [4]. The method has been 
widely used to simulate various problems in the complex fluids. For example, simulation 
of solute-solvent dynamics [1, 9], claylike colloids [5], star polymers in solution [11], and 
flow-induced polymer translocation [10]. 

As far as we know, MPCD method has never been used in the atmospheric hydro¬ 
dynamics or hydrostatics problems. So, the method will be tested for its ability when 
applied to atmospheric problems. For this purpose, MPCD will be applied to a simple 
problem, simulation of behavior of atmospheric neutral gas in hydrostatic equilibrium. 
The standard of hydrostatic atmospheric model used in this study will be described in 
the section 2, while details of the method and simulation mechanism will be explained 
in section 3. 


2 Standard Atmospheric Model 

Based on U.S. Standard Atmosphere 1976 [3, pp. 6-7], the atmosphere up to altitude 
86 km is assumed to be dry, homogeneously mixed, and obeys ideal gas law. The effect 
of homogeneity leads to composition of the atmosphere up to this altitude have uniform 
molecular weight, Mq = 28.96443 kg kmoR^. The atmospheric composition is dominated 
by nitrogen gas N 2 with 78.084 %, oxygen gas (20.9476 %), and traces of Ar, CO 2 , Ne, 
He, Kr, Xe, CH^, and H 2 . 

Atmospheric hydrostatic state can be described by using ideal gas law. Supposed 
atmospheric temperature T, particle density p, and molecular weight Mq, the atmospheric 
pressure is given by 


pi?*r 

Mo * 


( 1 ) 


Profile of particle density for any altitude can be deduced from equation (1) and is given 
by 


p{z) = Po exp 


H F 


( 2 ) 
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where zq, z are initial and current altitudes, respectively. Scale height ff is defined by 
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where i?*, g are universal gas constant and gravitational acceleration respectively. If we 
assume temperature T = Tq = 288.15 K, then the scale height is 8.4334 km. 


3 Multi-particle Collision Dynamics 


Simulation step in MPCD algorithm is divided into two steps, streaming and collision. 
In the streaming step, particle position is updated based on its velocity calculated in the 
collision step. The particle position for any time step is calculated by using 

fi{t + St) = fi{t) + Vi{t)St + 

In this equation, i express the index of the particles, f{t + 6t) is position at time t + 6t^ 
while r(t), v{t) are position and velocity at time t respectively. 

In the collision step, particles are sorted into a collision cell with lateral size a. MPDC 
algorithm allows the particles to interact with each other in the same cell only. The 
interaction itself is represented by stochastic rotation matrix with random angle for 
every time step. Particle velocity after interaction is calculated by using 

Vi{t + 6t) — u + R5f7, 6v = u — Vi{t)^ 


where f;, fl, R, 6v are velocity after interaction, velocity of center of mass, rotation matrix, 
and relative velocity between particle and center of mass, respectively. The velocity of 
center of mass is expressed by 
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where Nc is number of particle per cell. Two-dimensional rotation matrix is calculated 
by using 


R = 


cos{n6) — sm{n0) 
sm{n9) cos{n9) 


where n is chosen between +1 and —1 randomly, while 9 in this study is set to 90°. 

Interaction of the particles with the box is applied through the implementation of 
boundary conditions (BC). Periodic BC is applied on the left and the right of the box, 
while specular reflection and bounce-back is applied on the top and the bottom side 
respectively. Implementation of the last two BC is implemented by modifying tangential 
and normal components of particle velocity based on 


vt = (2r -1) vt, 

V = Vt + Vn, 
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where F = 1 for specular reflection and F = 0 for bounce-back. 

In the equation (2), we can see that the density profile is dependent only on altitude. 
Hence, motion of the particles in two-dimensional coordinate is sufficient to achieve this 
profile. In this study, The simulation is performed by assuming that there are N particles 
with mass m = 1, position r and velocity v. The particles are distributed randomly in 
two-dimensional box with size 10a x 20a, a = 1. The number of collision cells in the box 
is 200 and average number of particles per cell is fixed to 5. We set i?*T = 1, this is led 
to mean particle speed v = 1.5958, and speed of sound Cg = 1.1832. 

Each particle is given a random initial velocity which follows normal distribution. The 
distribution is generated from Box-Miiller transformation and expressed by [7] 

yi = A/-21og(Mi)cos(27rM2), 
y 2 = a/- 2 log(iti) sin(27rit2), 

where i/i, U 2 are uniformly distributed random numbers, and yi, y 2 are a pair of normal 
distributed random numbers. The pair of values are used interchangeably to produce the 
initial velocity for each particle through 


Vi^ya + 11] 


J ^ 1 , i even 
\ y 2 , i odd ’ 


where a and fi are standard deviation and mean of the distribution. In this study, /i = 0, 
while a is the thermal speed of particle. 

Mean free path is set to 0.25a and corresponding to collision interval 


6t = Ato5 


to = a 


Mo 

i?*r’ 


where to is unit time. Mean free path is taken smaller than lateral size (A <C a) indicates 
that there are no molecular chaos and Galilean invariant [6]. This problem can be solved 
by using grid shifting method [6]. The method is implemented by multiplying the grid 
position before collision with any random number b and then moving back after collision. 
The random number is chosen in the range +<^/2 and -<^/2. 

Since the system depends on value of gravitation acceleration, it will be varied in 
performing the simulation. The value is chosen among 1, 1.5, 2, 2.5, 3, and 3.5. The 
simulation itself is carried out until 5000 time steps. The density profile is generated from 
average density value from 2500 to 5000 time steps, while the particle speed distribution 
is taken from the last time step. 


4 Results and Discussion 

Simulation program is built by using Xojo (http://www.xojo.com) which based on BA¬ 
SIC language. The application is programmed and executed in Arch Linux with processor 
Intel i5 @2700 MHz and 10 GB DDR3 RAM. Figure la shows the average of simula¬ 
tion time for several total number of particles during both streaming and collision steps 
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for each time step. Although, We can not compare this result with simulation time of 
others methods, MPCD is reliable as an alternative of lightweight method for fluid sim¬ 
ulation. On the flgure lb. We can see that the algorithm satisfles energy conservation 
and therefore also fulflll linear momentum conservation. 



Figure 1: (a) Average simulation time for each total number of particles; (b) Profile of 
kinetic energy during simulation: g — 1 (dashed), g = 1.5 (solid), g — 2 
(dotted), g = 2.5 (dash-dotted), ^ = 3 (t), and g — 3.5 (♦). 


Figure 2 shows that the system has achieved equilibrium when time step is close to 
2000. This equilibrium is indicated by the trend of Mean Absolute Error (MAE) between 
density profile calculated from the model (2) and the simulation result. Figure 2 also 






Figure 2: Mean Absolute Error (MAE). 

shows that for small gravitational acceleration value {g < 1.5), fluctuations still occurred 
although the system has achieved equilibrium. Especially for g — 1 where the system is 
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not stable enough although the trend is quite stable. In contrast to that, more stable 
fluctuations are shown by higher values of g. 

According to Figure 3, MPCD algorithm is fairly accurate in modeling hydrostatic 
atmospheric system. The dash-dotted line on the flgure is calculated from equation (2), 
while the other one is the simulation result. Hydrostatic density profile can be achieved 






(d) g = 2.5 (e) ^ = 3 (f) g = 3.5 


Figure 3: Average density profile as a function of altitude after equilibrium: Theoretical 
prediction (♦ dot-dashed line) and MPCD results (• solid line). 


when g >2. However, when g > 3.5 it seems that the applied force is too big and attracts 
the majority of the particles to the bottom of the box, while for smaller particles 
are not really affected by applied force. These effects are also shown in the profile of 
particle speed distribution (Figure 4). Particle speeds do not obey Maxwell-Boltzmann 
distribution when gravitational acceleration is higher than 2.5. 


5 Conclusions 

MPCD algorithm is quite fast and can be used as an alternative of lightweight (low per¬ 
formance computing) method for fluids simulation. Atmospheric hydrostatic equilibrium 
can be achieved by using gravitational acceleration that is proportional to the particle 
mass. Based on the results, the recommended gravitational acceleration value is about 
1.5 < ^ < 2.5. MPCD method is capable to modeling atmospheric hydrostatic state and 
probably can be extended to hydrodynamic cases in atmospheric layer. 
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(d) g = 2.5 (e) 5 = 3 (f) g = 3.5 

Figure 4: Particle speed distribution: MPCD result (bar) and fitted (solid line). 
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